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Abstract 

Streamers are the first stage of sparks and lightning; they grow due to a strongly enhanced 
electric field at their tips; this field is created by a thin curved space charge layer. These multiple 
scales are already challenging when the electrons are approximated by densities. However, elec- 
tron density fluctuations in the leading edge of the front and non-thermal stretched tails of the 
electron energy distribution (as a cause of X-ray emissions) require a particle model to follow 
the electron motion. As super-particle methods create wrong statistics and numerical artifacts, 
modeling the individual electron dynamics in streamers is limited to early stages where the total 
electron number still is limited. 

The method of choice is a hybrid computation in space where individual electrons are fol- 
lowed in the region of high electric field and low density while the bulk of the electrons is 
approximated by densities (or fluids). We here develop the hybrid coupling for planar fronts. 
First, to obtain a consistent flux at the interface between particle and fluid model in the hybrid 
computation, the widely used classical fluid model is replaced by an extended fluid model. Then 
the coupling algorithm and the numerical implementation of the spatially hybrid model are pre- 
sented in detail, in particular, the position of the model interface and the construction of the 
buffer region. The method carries generic features of pulled fronts that can be applied to similar 
problems like large deviations in the leading edge of population fronts etc. 
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1. Introduction 



1.1. X-ray bursts and terrestrial gamma-ray flashes 

Terrestrial Gamma-Ray Flashes were first observed accidentally in 1994 f^; later they were 
found to be correlated with thunderstorm activity. In 2001, energetic radiation was also observed 
from normal cloud-to-ground lightning and later from rocket triggered lightning jstl- Mean- 
while X-ray bursts were also found in the laboratory near long sparks flSSBl- Such flashes 
and bursts are likely to be produced by Bremsstrahlung of very energetic electrons (so called 
run-away electrons). 

One possible source of highly energetic electrons are the streamer discharges that pave the 
way of sparks and lightning. Streamers are ionized plasma channels that grow into a non-ionized 
medium due to the self-enhancement of the electric field at their tips. In this high field region, 
the electron energy distribution is very non-thermal and can have a long tail at high energies; 
it could act so much as a self-organized local electron accelerator that it is a possible source of 
X-rays 

The single electron dynamics in streamers can be studied by a particle model, which models 
the streamer dynamics at the lowest molecular level. It follows the free flight of single electrons 
in the local field between abundant neutral molecules and includes the elastic, exciting and ion- 
izing collisions of electrons with the molecules in a stochastic Monte Carlo procedure. However, 
the increasing number of electrons in a growing streamer channel leads to an increasing demand 
of computational power and storage and excludes the simulation of all electrons in a full long 
streamer on present day's workstations. 

The generation of highly energetic run-away electrons from streamers therefore has only been 
modeled in particle models with simplifying assumptions or tricks. In |8;] a Monte Carlo simu- 
lation treated planar fronts, and 3D aspects were modeled by a simplified electric field profile in 
the forward direction; in llOll a 3D axis-symmetrical streamer was simulated with super-particles 
of high weight where many real electrons are replaced by one super-particle. Both simulations 
indicate that electrons can be accelerated to very high energies in the high field zone at the tip of 
streamers. However, the ID simulation |8|1 does not include the intricate 3D multiscale spatial 
structure of the streamer head properly ilUll2lll3[l . while the super-particle approach jlOtl has 



been shown |14l] to lead rapidly to numerical artifacts already shortly after a streamer emerges 
from an ionization avalanche; it will not represent the physically important tail of the electron 
energy distribution correctly. 

1.2. Need and concept of spatially hybrid computations for streamers 

Realizing that the number of electrons in the high field zone at the streamer tip is limited, 
and that only these electrons have a chance to gain high energies, a natural idea is to include 
only those electrons into a particle model while all others are treated in density or fluid approx- 
imation, as the fluid approximation is cornputationally much more eflicient and has been widely 
used in streamer modeling, see e.g., ful |H fisj [l6l |H fll fl9i |20i [III 



Such a spatially hybrid approach should allow us to follow single electrons in the relevant re- 
gion while avoiding the introduction of super-particles. We here present essential steps of its 
implementation. 

The concept of a spatial coupling of density and fluid model is shown in Fig.[T] The natural 
structure of the streamer consists of an ionized channel and an ionization front. In the chan- 
nel, the electrons are dense enough to be approximated as continuous densities, and the field is 
too low to accelerate them to high energy. In the ionization front the electron density decreases 
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Figure 1 : (Color online) The spatially hybrid scheme is shown in a one dimensional section of the streamer The solid 
line shows the electron density ng as a function of the spatial coordinate z- The particle model follows the leading part of 
the ionization front and the fluid model is applied to the rest. The questions treated in the present paper are: Which fluid 
model approximates the particle model in the best manner? Where should the model interface be placed and how should 
it be structured? The figure is reproduced from j^l ■ 



rapidly while the field increases, and the electrons gain high energies far from equilibrium. The 
hybrid model should apply the particle model in the most dynamic and exotic region with rela- 
tively few electrons and the fluid model in the remaining region that contains the vast majority 
of the electrons. The hybrid model is not only suitable for studying possible run-away electrons 
from streamer heads, but also for studying the role of electron density fluctuations on streamer 
branching, and for studying the gas chemistry inside the streamer and the electron energy distri- 
bution at the streamer tip in general. 

To implement such a hybrid computation, we have set first steps in |@, 28 ] . In |@] we studied 



the Monte Carlo discharge model in detail and compared it with its fluid approximation; we did 
that for electron avalanches in a constant electric field and for planar ionization fronts; the study 
was for negative streamers in nitrogen. We derived reaction and transport coefficients for the 
fluid model from Monte Carlo electron avalanches, and we found that the mean electron energy 
increases at the front of an electron swarm, even if the electric field is constant. As the local 
field approximation cannot break down in a constant field, we attributed this effect to a break- 
down of the local density approximation of the fluid model in the steep density gradients of the 
advancing front. In the fast track comunication 12811 . we presented a first realization of a spatial 
coupling of particle and fluid model for planar ionization fronts. However, a closer investigation 
of these results indicates that further improvements are desirable, and the paper did not include 
details of our work. Therefore we here present details and further investigations on how to couple 
particle and fluid model in a setting of planar fronts. This lays the basis for full three dimensional 
streamer computations in the near future where the dynamics of single electrons in the relevant 
region can then be fully followed. 



1.3. Particle fluctuations and large deviations in pulled fronts 

The hybrid coupling in space that we present here, is designed for the simulation of streamer 
discharges. However, the method carries generic features of so-called "pulled fronts". The 
pulling terminology was probably first used in 1976 in the biomathematical literature ll29ll for a 
phenomenon found in 1937 independently by Kolmogoroff, Petrovsky and PicounofF and by R.A. 
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Fisher, examples include chemical bacterial |31], and virus invasion fronts |32]. "Pulling 



means that the penetration of one phase or species into another is dominated by the instability of 
the penetrated state; front shape and velocity are completely determined by the linear instabil- 
ity of the penetrated state, and not by a balance of forces of two competing states (for reviews 



see Il33ll34ll '). Therefore, if the penetrating species is organized in particles (be it plants, bacteria, 
molecules or electrons), the penetration of one particle into the unstable state can create a local 
avalanche, and it can largely influence the dynamics. Large deviations from the mean in the 
leading edge of a front can therefore grow out into macroscopic structures. Examples are the 
first electron that reaches a completely non-ionized region with a local field above the ionization 
threshold, or the first species that reaches a new habitat 1 35 , 3^ 37 ] . In that sense, it is surprising 



how well fluid approximations work for many pulled fronts. 

Spatially hybrid models therefore are of general interest and a viable method for pulled front 
problems as they treat the dynamically relevant leading edge with low particle densities in a 
particle model and the bulk of particles behind the front in density approximation. In particular, 
for off-lattice models it poses a challenge. Questions include: Which fluid model approximates 
the particle model in the best manner? Where should the model interface be placed and how 
should it be structured? 

Adapting the buffer zone method H 111 lil, H to connect the particle with the fluid 
zone, we encounter two more features that are generic for pulled fronts. 1. When the model 
interface moves along with the velocity of the front, the particles in this frame of reference move 
backward. This is because the front is pulled by motion and multiplication of particles in the 
leading edge where the electric field and the electron drift velocity are the highest, therefore 
throughout the front, individual particles on average move backward, even in the leading edge. 
As we will show in Section 13.41 this leads to a large technical advantage for the construction 
of the buffer region, as particles do not need to be created out of the fluid region, but they only 
have to be absorbed there. 2. Avalanches in the unstable region have essentially the same front 
shape and velocity as full fronts 1331 19I1. When particle avalanches are used to derive transport 
and reaction functions for the fluid approximation, a pulled front is therefore well matched by 
this model. 

1.4. Content of the paper 



Section |2] discusses why our first hybrid coupling approach 12811 still showed discontinuities 
at the model interface. The answer is found in inconsistent fluxes, and the problem is resolved 
by improving the fluid model through a gradient expansion in the density. Then transport and 
reaction coefficients are derived for this extended fluid model, and the extended fluid model is 
compared to the full particle model, both for avalanches and for planar fronts. The extended fluid 
model now describes the dynamics well, except for particle fluctuations in the leading edge of the 
front. The coupling of the extended fluid model with the particle model is described in Section|3] 
where the coupling concept and implementation are given in detail; discussed are the hybrid 
algorithm, the numerical implementation, the position of the model interface and the structure 
of the buffer region at the interface. While this whole study is done in a fixed background field 
for negative streamer fronts in nitrogen, Section|4]investigates how these results are generalized 
to other fields. Section |5] contains the conclusions. Appendix lAl discusses an alternative to the 
extension of the fluid model by a gradient expansion, namely the extension with an electron 
energy equation. 
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2. Extending the fluid approximation 



The particle model follows the individual electrons on their deterministic flight in the lo- 
cal field until the next collision with a neutral particle. Time and kind of collision and energy 
and momentum of the electron after the collision are calculated in a Monte Carlo procedure, in 



accordance with the cross-section data base of Il44ll . Details of the procedure are explained in 



The classical fluid model as used by many authors for negative streamers lU, 12, 15, 3 12 



can be derived from the particle dynamics through the Boltzmann equation 
assuming the local field and the local density approximation, where the electron mean energies, 
transport coefficients and reaction rates are functions of the local reduced electric field and the 
local electron density. (We remark that the local field approximation is typically emphasized in 
this derivation, while the local density approximation is not.) However, we have shown that in 
electron avalanches, the local mean energies of the electrons vary from the tail to the front of the 
swarm Jgt], even though the electric field is constant. It is impossible to improve the fluid model 
to the same level of description as the particle model, and we do not intend to do that. But the 
main features that create inconsistencies and that influence the coupling of the two models, have 
to be investigated and understood. 

Here we investigate why the electron fluxes are inconsistent when coupling the classical fluid 
model with the particle model, and we extend the fluid model with a gradient expansion in the 
electron density to incorporate nonlocal effects. The two fluid models are compared with each 
other for both the electron avalanche and the planar front. 

2.1. A short recollection of the classical fluid model 

To elucidate the essential structure of the problem, we here analyze discharges in pure nitro- 
gen. The classical fluid model then contains two continuity equations for the densities yig, iip of 
electrons and positive ions 

^ + V.j, = Si, (1) 



dt 



Si, (2) 
-M£)E«.-D,(E)-V«„ (3) 



where is the electron flux and S is the source of electrons and ions, /i; represents the mobility 
and D; is the diffusion matrix, E and E are the electric field and the electric field strength; the sub- 
script "/" stands for "local" and denotes the classical fluid model in contrast to the extended one 
to be introduced below. The densities of charged particles together with the boundary conditions 
determine the electric field through the Poisson equation 

e («„ - lie) 

V.E = — -. (4) 

The source term accounts for the impact ionization in local field and local density approximation 

Si = |«, ME) E| aiiE). (5) 
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Figure 2; Hybrid simulation of a planar front propagating in = -100 kV/cm with the model interface at -60 kV/cm. 
The hybrid model introduced in \\% is used. Electron density ("+") and electric field (dot-dashed line) are plotted, and 
a discontinuity in the electron density is marked with a dashed circle. 



2.2. Why the fluid model needs an extension 

2.2.1. Discontinuities in the hybrid model for planar fronts 

We have derived transport and reaction coefficients for the classical fluid model from the 
particle model in |0], and we have shown first results of a hybrid computation which couples the 
classical fluid model with the particle model in lEsll . 

However, to ensure a stable and correct interaction between the two models, we later inves- 
tigated the electron flux densities on the model interface, and the mean electron energies and 
velocities around the model interface, and we found a disagreement of the local electron flux 
between the two models. When the particle and the classical fluid model are applied in the same 
region of the ionization front, the mean electron flux density is lower in the particle simulation 
than in the classical fluid simulation as we will demonstrate below. If one places the model 
interface near the maximum of the electron density, a density jump can appear near the model 
interface. This is shown in Fig. |2] where the model interface is placed at the electric field E - -60 
kV/cm with the field ahead of the front being - -100 kV/cm; here the local fluxes are dis- 
continuous across the model interface, and this results in the visible discontinuity of the density. 
The numerical discretization and other details of this experiment can be found in 128 1. 

The effect appears although the transport and reaction coefficients were generated in the par- 
ticle swarm experiments and then implemented into the fluid model; therefore the two models 
should be consistent. To understand the underlying problem, we review the conceptual differ- 
ences between the two models and compare them in a electron avalanche experiment. 



2.2.2. Electron avalanche simulations in particle and fluid model - a reinvestigation 

Fig.|3]shows simulated electron swarms propagating in a uniform constant field of -lOOkV/cm. 
The simulations have been carried out both with the particle model (solid lines) and with the clas- 
sical fluid model (dotted lines) in ID starting from the same initial conditions. For the particle 
model, ID means that the simulations have been carried out in a volume with a narrow transversal 
cross-section, with periodic boundary conditions at the lateral edges. 
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Figure 3: Evolution of electron swarms at three consecutive time steps ti, ti, f3 in a field of —100 kV/cm, where fi , fi, 
and f3 are 20, 130, 240 ps when the simulation starts from 100 pairs of electrons and ions. Shown are the spatial profiles 
of the total swarm in the ID particle model (solid Hne) and in the ID classical fluid model (dotted line). Also shown are 
the spatial profiles of those electrons that were already present at time 1i (at the same three time steps). 

Here and in all later comparisons of particle or fluid results, the initial distribution of electrons 
and ions is generated in the following manner An electrically neutral group of 100 electrons and 
ions is inserted at one point in space. Their temporal evolution is calculated with the particle 
model. After a short time that depends on the electric field, e.g., after ti - 20 ps in a field of 
100 kV/cm, a small swarm forms in which the electron density profile is well approximated by 
a Gaussian distribution. This distribution of electrons and ions is used as an initial condition for 
all simulations, including particle, fluid and hybrid. 

In the particle model, the standard Monte Carlo technique in a constant electric field is ap- 
plied. For the fluid model, the equations of the particle densities are discretized in space with 
a finite volume method. The particle densities are updated in time using a third order upwind- 
biased advection scheme combined with a two-stage Runge-Kutta method. The same time step 
and cell size were used in the particle and the classical fluid calculation, i.e., Af = 0.3 ps and 
Az = 2.3 um. More details of the numerical implementation of particle and fluid model can be 
found in |9fl . 

In the avalanche experiment, we let both the particle model and the classical fluid model 
follow the swarm from t\ — 20 ps to f3 - 240 ps, and the electrons in the particle simulation 
are mapped to densities on the same grid as used in the fluid simulation. During the simulation, 
we follow the growth of the total electron avalanche, and we pay specially attention to the group 
of electrons that are initially (at time ti) present as seed electrons. Fig. [3] shows the electron 
density profiles of the swarm at time t\ and at the later stages t2 and fa, and also the density 
distributions of only those electrons that already existed at time t\, at all three time steps ti, t2 
and f3, neglecting all electrons generated later 

Fig. [3] shows that the profiles of the total swarm are nearly the same in the classical fluid 
model and in the particle model. This was to be expected as we just derived the transport and 
reaction coefficients ///, D/, and or/ for the fluid model through this condition |j9j[|. However, 
following only the swarm of electrons that already existed at time t\ , the particle swarm stays 
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behind the fluid swarm. This observation has to be interpreted in the foUowing manner In the 
fluid simulation, the electron swarm grows homogeneously, i.e., the center of mass of the initially 
present electrons moves with the same speed as the center of the mass of total electron swarm. 
However, in the particle simulation the mean displacement of the whole swarm is larger than that 
of the initially present particles. This effect can only be caused by a larger growth rate of the 
electron swarm at the tip than at the back of the swarm. This shows that the classical fluid model 
is based on approximations that here become inaccurate. 

The electron flux rates in particle or fluid model can be derived by averaging the local electron 
flux je over the simulation domain of volume V and over a short time interval t 



where A^^, is the total number of electrons. We find a higher mean flux density per electron in 
the fluid than in the particle swarm simulation. The mean flux density per electron equals the 
mean electron drift fiE, as the difi'usive flux vanishes if all electrons are inside the integration 
volume, see Eq. (O. The higher electron mobility in the fluid model is consistent with the larger 
displacement of the swarm of initially present electrons. 

We conclude that the particle fluxes are inconsistent between particle model and classical 
fluid approximation both for planar fronts and for electron avalanches. For the avalanches, this 
has been shown both by interpretation of Fig.[3]and by a calculation of fluxes through Eq. (|6|. 

2.2.3. Discussion: validity of the local density approximation 

That both the mean flux density Q and the avalanche growth characteristics (Fig.O differ 
between particle model and classical fluid approximation, is actually due to the same reason: the 
local density approximation. As the electric field is constant in the swarm experiment, it cannot 
be the local field approximation. 

In 101, it was already shown that the local electron energy distribution depends on the local 
electron density gradients both in planar fronts and in electron avalanche. It indicates that even 
in a uniform constant field, the ionization rates increase from the tail to the head of the electron 
swarm. Therefore there are two contributions to the speed of the swarm: i) the electron drift in 
the local field as a main contribution and ii) the unequal growth of the swarm due to the unequal 
electron energy distribution. 

The classical fluid model assumes that the local mean energies and local mean reaction rates 
are functions of the reduced electric field, and it ignores the fact that within the same electric 
field, the ionization rate can be different. At the front of the swarm, this leads to a lower reaction 
rate in the fluid simulation than in the particle model. But by overestimating the ionization rate 
at the tail, the swarm in the classical fluid model generates the same amount of electrons as in 
the particle swarm. And by overestimating the mobility of electrons, the swarm nevertheless 
propagates with the same speed. Or in other words, the classical fluid model approximates the 
swarm density profile rather well, but at the price of wrong local fluxes and wrong local reaction 
rates. In simulations with non-uniform fields, such as streamers in which the density profiles are 
determined only by the front while the tail doesn't play a role in the propagation, the classical 
fluid model pays the price by generating lower electron densities inside the streamer than the 
particle model |@] . 

Therefore the fluid model has to be extended to include the dependence of the ionization rate 
on non-uniform electron density conditions. And by introducing a nonlocal term, we expect the 
fluid model to supply the same electron flux as the particle model within the same field. 




(6) 
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2.3. Gradient expansion and extended fluid model 
2.3.1. The introduction of the extended fluid model 

The electron ionization rates for non-uniform fields and electron densities are discussed 
in 13 47]. There perturbation theory is used to obtain the rate coefficients, and it is shown 



that they can be represented by a gradient expansion about the local values 

= .S, (l + /tie ■ V In + /tzA^/EV ■ (E/A^) + ■ V ln(£/Af)), (7) 

where n^ is the electron density, E and E are the field and the field strength and e = E/Zs is the 
unit vector in the direction of the electric field, is the molecule density of the background gas, 
and liQ, k\, k2, k^, are parameters depending on E/N that have to be determined. 

In 10], we have shown that the local density approximation is insufficient while the local field 
approximation is not problematic. Therefore, we neglect the field gradient terms, focus on the 
density gradient in Eq. and approximate the source term in our nonlocal model for fixed 
density as 

S„ = S,[l+ki{E)e-W\nn,) 

= IJ„{E)a„(E){En, + ki(E)E-Vn,), (8) 

where the subscript n denotes the nonlocal or extended fluid model. Otherwise the model is 
identical to the classical fluid model, except that all transport and reaction coefficients have to be 
determined a new, they are denoted as jj„, D„, and a„. All these functions as well as ki are now 
determined. 

It has been known for long that the consistency of the fluid model with particle swarm results 
is important. Robson et al. 1480 recently reviewed this issue and emphasized that the consistency 
requirement actually applies to transport coefficients and reaction rates. In l^f], we therefore have 
determined the mobility fiiiE), diffusion tensor D;(£'), and ionization rate ai{E) for the classical 
fluid model from particle swarm experiments. 

2.3.2. Electron mobility fi„ 

When the classical fluid model is extended with a density gradient term, the transport and 
ionization coefficients may also change and should be re-defined from the swarm experiments. 
As discussed above, the mobilities can be determined either i) from the mean displacement of 
electron swarm in a uniform constant field 

= -, 

h - h 

(as was done in the classical fluid model) or ii) from the mean displacement of the initially 
present particles 

z'-, -zi 
h - 1\ 

where z\ and z-i are the centers of mass of the swarms at times t\ and f3, and z'-^ is the center of 
mass of those electrons at time t^ that were already present at time t\, cf. Fig. [3] The second 
definition represents the real electron motion and is used in the extended fluid model; within 
numerical accuracy, this mobility coincides with the one derived from the average local flux rate 
On the contrary, the first definition is actually the sum of the local electron fluxes under the 
influence of the electric field plus the nonlocal growth of the swarm. 
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Figure 4: fii and )i„ as functions of the electric field strength E. 



The mobilities /i/ or /i„ are plotted in Fig.H] Here jU/ is the electron mobility in the classical 
fluid model and ju„ is the one in the extended fluid model. They are almost the same for fields 
E below 30 kV/cm when impact ionization is small or even completely negligible, but the dif- 
ference increases as the field strength increases. It can be noted that the curves in Fig.|4]are not 
completely smooth. This is due to stochastic fluctuations within the Monte Carlo simulations. 
However, it has been checked that these small fluctuations have a negligible influence on the final 
results. 



2.3.3. Impact ionization rate a„ 

The reaction rate a is derived from the growth rate of the total electron number Neit) in a 
particle swarm simulation. Within a time r, the swarm propagates a distance ii{E)Et, and the 
reaction rate is determined by 

lnjV,(f2)-lnjV,(fi) 
H{E)E (t2 - ti) 

When is replaced by /^„, the ionization rate changes from a; to 

a„{E) = a,{E) (10) 
as the simulation fixes the product jjiai - ii„a„. 

2.3.4. Diffusion tensor D„ 

The diffusion rates can also be obtained according to two different definitions: the diff'usion 
of the electron swarms D/ or diffusion of the initially present electrons D„. The results turn out 
to be the same 

D,(E) = D„(E), (11) 

as shown below. 
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Figure 5 : The parameter function k[ as a function of the electric field. 



2.3.5. The parameter function ki of the gradient expansion 

The extended fluid model needs die function ki in Eq. dS). The rates ^i, k2, and k^ in Eq. d?) 
were calculated by Aleksandrov and Kochetov |!46ll by solving the Boltzmann equation in a two- 
term approximation. But to stay consistent with the approach above, all parameter functions in 
the fluid model are here derived from particle swarm simulations. Since we neglect the field 
gradient term in the source term, only ki needs to be calculated. 

In fact, ki can be calculated by simply comparing the electron continuum equations in both 
the classical and the extended fluid model. As both fluid models can properly describe the swarm 
when appropriate parameters are applied, we have 

d,ne = V ■(///(£■)£«,) + V ■ (d,(E) ■ V«,) + «,///(£)£■ ff/CE) (12) 
= V ■ + V ■ (d„(E) ■ V«,) + n, E a„(£) + ki{E) p„{E) a„(E) E ■ V«,. 

In a constant uniform electric field such as in the electron avalanche experiment, lUi(E), fi„(E) 
and E are constant and V ■ (/i/(E) E «p) = /i/(E) E ■ Vn^, etc. Removing identical terms on both 
sides of the equation, we retrieve D/ = D„ ( fTTT i and we get fii(E) - ju„(£')(l + ki{E)an{E)) or 



ki{E) 



Idi(E) - fi„(E) 
a„(E) fi„{E) 



(13) 



where /ii, ju„ and a„ were derived above from the particle swarm experiments, as a function of 
the field is shown in Fig|5] 

2.4. Comparison of the extended fluid model with the particle model 

Now the stage is set to compare the extended fluid model with the particle model both in 
swarms and in planar ionization fronts. The extended fluid equations are discretized in the same 
manner as the classical fluid model, and the particle densities are updated with the same scheme 
as before. The same time step and cell size are used in the extended fluid calculation as in the 
particle and the classical fluid calculation (see Sect. 12.21 ). 
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Figure 6: The density profiles of electron swarms and the initial electrons in a field of -100 kV/cm are shown at times 
ti , t2 and . The plot is the same as in Fig.|3] but now the extended fluid simulation is added as a dashed line. 

2.4.1. Swarm simulations 

In Fig|6] we show electron swarms at times ti, t2, and f j in a constant field of -lOOkV/cm; the 
swarms were followed by particle simulation (solid line), classical fluid simulation (dotted line) 
and extended fluid simulation (dashed line). For the whole swarm, all three models give similar 
results, but the extended fluid model follows the evolution of the initially present electrons much 
better than the classical fluid model. 

However, the figure also shows that in the leading edge of the swarm (marked with a large 
dashed circle) neither the extended fluid model nor the classical fluid model describe the electron 
density distribution of the particle model precisely. As described in jstl, individual electrons with 
high energies largely deviate from the mean, and none of the two macroscopic fluid models can 
reproduce this microscopic behavior. The densities in this region are 3 to 5 orders lower than the 
maximal density. But as streamer ionization fronts are so-called "pulled" fronts, cf. section [131 
the behavior in the leading edge of the front actually determines the front velocity. From this 
point of view, the hybrid model that uses the particle model in the leading edge of the ionization 
front, is the method of choice. 

2.4.2. Front simulations 

Fig. |7]shows the temporal evolution of the planar front in a field of = -100 kV/cm in 
the particle simulation (solid line), classical fluid simulation (dotted line), and extended fluid 
simulation (dashed line). Compared to the classical fluid model, where the maximal electron 
density in the front and the saturation level of the ionization behind the front are about 20% 
lower than in the particle model, the extended fluid model approximates the particle model much 
better. The particle and fluid front move with approximately the same velocity, but the particle 
front moves slightly faster than the extended fluid front, and the extended fluid front moves 
slightly faster than the classical fluid front, in agreement with Fig.|6l (We recall from [21 that the 
leading edge of a swarm and of a pulled front in the same field have the same spatial profile and 
create the same velocity.) 
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Figure 7: Temporal evolution of the electron densities in a planar front in a field of E'^ = -100 kV/cm. Shown are 
the spatial profiles of electron densities derived with the particle, classical fluid or the extended fluid model at time 
levels t=0.672 ns, 0.816 ns and 0.96 ns (particle simulation: solid line, fluid simulation: dotted line, and extended fluid 
simulation: dashed line). 




Figure 8: The relative density difl'erence )14t between extended fluid model and particle model. 
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Having analyzed planar fronts at - 100 kV/cm, we now summarize the front results for fields 
ranging from -50 to -200 kV/cm in Fig. [8] The figure shows the relative difference 

^e,part~^e,e.fluid 

1 — • (14) 

tl 

e.part 

of the saturated electron density behind the front in the particle model (n^^^,.,) and the ex- 
tended fluid model («^^^ /^/,„v/)- Compared to the relative difference between the particle model 
and the classical fluid model discussed in which increases from 10% at 50 kV/cm to 40% 
at 200 kV/cm, we now find that the relative density diflTerence in the extended fluid model never 
exceeds 10% within this range of fields. 

2.5. The drawbacks of the extended fluid model 

By approximating the nonlocal ionization rate with a density gradient expansion, the ex- 
tended fluid model reproduces the ionization level behind the front much better than the classical 
fluid model. However, Fig.[8]shows that this ionization level in the extended fluid model exceeds 
the one in the particle model when the field E is below 125 kV/cm. This will not harm its cou- 
pUng with the particle model, but the reason for this unexpected behavior is briefly discussed 
here, because it could lead to further improvements of the fluid model. 

Several possible reasons have been examined, for example, the quality of the parameters 
(ju, D, a and ki) used in the extended fluid model, the discretization of the gradient expansion 
term, and the relative ionization rate in the very leading edge where the electron density simply 
vanishes in the particle model, but remains nonzero in the fluid model due to electron diffusion. 
The answer was found by comparing the ionization rates within the front between both models. 

The simulation is done as follows. We first let the particle model follow the evolution of an 
ionization front until it is fairly smooth. This creates an initial condition that now is run further 
both with the extended fluid model and with the particle model. The two models are followed 
for a time of 100 Af where Af = 0.3 ps in a field of - -100 kV/cm. The result is shown in 
Fig- 121 In the first panel we plot the electron density profile in steps of 20 Ar in the particle model 
(dashed) and in the extended fluid model (solid). Of course, the curves are identical initially, 
and a difference builds up in time in the high density region of the front (marked with a dashed 
circle). In the second panel, we show the electric field which does not vary much between the 
two models during 100 Af. In the third panel, we show the relative growth of the ion density 
(drnp)/ne integrated over a time of 20 Af, which is identical to the local ionization rate S/rie 
within this time according to Eq. (|2|. The figure shows that the ionization rate is slightly higher 
in the extended fluid model in the region of large rig (marked with dashed circle), and lower in 
the region of small Note that these differences appear while the electric fields are the same. 

The analysis shows that the extended fluid model does not completely reproduce the particle 
model. This is understandable. The actual ionization rates S/n^ depend on the energy distribution 
of the local electrons. The gradient expansion in the ionization rate relates this energy distribution 
to the local field and to the electron density gradient. The single adjustable function ki{E) in this 
gradient expansion is chosen in such a way that an electron swarm in this field is well fitted. 
But an electron swarm has a characteristic spatial profile in a given field. Regions inside the 
front might combine a given field with a different density gradient for which the model has not 
been adjusted. A solution would be to allow for more adjustable functions inside the fluid model 
by expanding in higher order gradients. In the end, only a fluid model with infinitely many 
adjustable functions would appropriately describe the averaged behavior of the particle model. 
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Figure 9: The ionization front propagates into -100 kV/cm in both particle model (dashed) and extend fluid model 
(solid). The two simulations start from the same front which is pre-produced by a particle simulation, and follow the 
front propagation for 100 At where At = 0.3 ps. We present the electron density (first panel), the electric field (second 
panel), and the ionization rate (third panel) in time steps of 20 At. 
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We therefore conclude that the density gradient expansion is a substantial improvement of the 
fluid model, but that it does not contain the full physics of the particle model. 

An alternative improvement of the fluid model next to a gradient expansion was suggested 
to us recently (see acknowledgement), namely the introduction of the energy equation. With the 
energy equation calculating the local mean electron energies, the fluid model can describe how 
the electron dynamics depends on the mean local electron energies rather than on the electric 
field 1 49, 13, 45, 51, 52, 53]. In appendix lAl the fluid model with energy equation is used 
to follow the electron avalanche in a constant electric field. Comparison of Fig. |6] and Fig. [15] 
shows that the extended fluid model follows the full particle avalanche a bit better, but that the 
fluid model with energy equation fits the mean electron energies in the particle avalanche quite 
well. 

We now proceed to the hybrid model where the extended fluid model will be used only in the 
less critical inner region of the streamer. 



3. The hybrid model 

Our goal is to build a fully 3D hybrid model for streamer channels, which couples the particle 
model with the extended fluid model in space as shown in Fig.[T] With the experiments carried 
out for planar fronts, we would like to test our hybrid concept, and find out how to apply the 
different models in suitable regions adaptively and how to realize a reasonable coupling of those 
models. Since most of the electrons will be approximated as densities in the hybrid simulation, 
the problem of the enormous number of electrons in a 3D particle streamer simulation can be 
reduced. Following the particles can still be a heavy burden for a hybrid model. Therefore 
the hybrid simulation will let the fluid model approximate most of the electrons as densities and 
leave as few electrons as possible for the particle model, under the constraint of producing similar 
results as the pure particle simulation. 

In this section, we present the algorithm for the ID hybrid calculation and discuss the numer- 
ical details. Two important issues of the ID coupling are discussed in detail: where particle and 
fluid model are to be applied, and how they should interact with each other 

3.1. The hybrid algorithm 

In our hybrid model, the particle model is used in the leading part of the ionization front and 
the fluid model in the rest of the domain. Between them, we have a model interface. The position 
of the model interface can be chosen either according to the electron densities or according to the 
electric field. This means that the model interface is either located at some electron density level 
He = X tie^max ahead of the electron density peak n^,,,,,, v, or at some electric field level E - y 
where E^ is the electric field ahead the front. In both cases, the real numbers x, y have to be 
chosen appropriately within the interval [0, 1]. We discuss this choice and the corresponding 
position of the model interface in detail in Sect. 13.31 

Suppose that the position of the model interface is located somewhere as shown in Fig. [1] 
As the streamer propagates forward, because the electric field moves with the front, the model 
interface moves with the ionization front because we keep it at - x nejnax or E - y E^ . In 
this way, the particle model at any moment of the simulation follows only a limited number of 
electrons in the low density region of the ionization front. Because the computational costs of the 
fluid model are small compared to the particle model, it is clear from Fig.[T]that we gain much 
efficiency in the hybrid model as compared to a pure particle simulation. 
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Figure 10: The flow chart of one hybrid computational step from to . 
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The actual gain due to the hybrid method will be even more pronounced in the case of fully 
three dimensional simulations. With this moving model interface algorithm, the hybrid model 
may simulate the full streamer dynamics in 3D without using super-particles while remaining 
computationally efficient. 

Once the position of the model interface is set, the interaction of the two models around 
the model interface needs to be established. To simulate the electron flux crossing the interface 
between the two models, we introduced a so called "buffer region" which is created by exten ding 
the particle region one or a few cells into the fluid region. It has been used in |3^ 4^ 41] 



for rarefied gases by coupling a Direct Simulation Monte Carlo (DSMC) scheme to the Navier- 



Stokes equation, and for other applications 11421 14311 as well. The buffer region helps to build a 



pure particle description around the model interface, such that the local particle flux in this region 
can be obtained as in a pure particle simulation. The particle flux across the model interface 
influences the total number of particles in the particle model, and it will lead to a corresponding 
increase or decrease of the density in the fluid model. 

In this way, around the area where the two models are coupled, the electron flux in a pure 
particle simulation is maintained in the hybrid computation, while global mass conservation 
holds. However, to obtain an accurate flux at this interface, the electrons in the buffer region 
near the interface should maintain a correct density profile and velocity distribution. In many 
cases 14111 the buffer region or a part of the buffer region at each time step have to be 



refilled or reconstructed with a large number of electrons with artificial distributions in energy 
and space. Such a filling or reconstruction on the one hand ensures a stable flux at the model 
interface, but on the other hand, it can create non-physical artifacts, which might cause wrong 
results. In our approach, the generation of electrons from a kinetic prediction can be avoided. 
The details are presented in Sect. 13.41 

Fig.[TO]shows the flow chart of the hybrid computations. At the beginning of each time step, 
the particle density is known in the fluid region and the individual particle information is known 
in the particle region and the buffer region. To update them to the next time step, we first move 
the particles in the particle and buffer region one step further The number of electrons crossing 
the model interface is recorded during the updating. This counted flux is used as a boundary 
condition to update the densities in the fluid region. To calculate the electric field, the particles 
in the particle region are mapped to the fluid grid as densities. In the end, we check whether the 
model interface has to be moved. A more detailed explanation is given in the following section. 

3.2. Numerical implementation 

In the particle model, the positions and the velocities of the electrons are updated with the 
leap-frog method. The electric field is solved on a uniform grid G (the fluid equations are dis- 
cretized on the same grid) with cells 



, / = 1,2,---M-, (15) 



where M- are the number of grid points with cell centers at z, = (/ - 1 /2)Az in the z direction, 
and Az is the cell size. 

The simulation begins with a few electron and ion pairs with a Gaussian density distribution. 
These initial particles are followed only by the particle model in the beginning of the simulation. 
As new electrons are generated, the number of particles eventually reaches a given threshold, 
after which the simulation switches to the hybrid approach. The threshold in our simulation 

18 



E (kV/cm) 50 



75 



100 



125 



150 



175 



200 



ly (jjm) 55.2 

L (mm) 6.9 

A^r (10'') 3 

Tt (ns) ~6.6 



41.4 
3.45 
3.2 
-1.2 



27.6 
2.76 
3.5 
-0.3 



23 
2.3 
3.8 
-0.18 



18.4 
1.84 
4 

-0.15 



13.8 
1.38 
5 

-0.12 



9.2 
1.15 
7 

-0.12 



Table 1: List of the transversal length /, (the transversal area is 1^ X ly), the longitudinal length of the system the 
threshold number Nt , and the time Tt for the particle simulation TO generatE Nr electrons for a number of electric 
fields. 

is normally set to be several million electrons to ensure a satisfactory statistics. With proper 
transversal area A of the system, the simulation with such number of electrons is already in the 
streamer stage, which means that a steady moving ionization front has developed and the charge 
layer at the streamer head totally screens the field inside the channel. In Table [T] we summarize 
the transversal areas (/,• x Ir) and the length of the system in front propagation direction L for 
varying fields, and we also list threshold numbers Nt and the time (Tt) for the particle simulation 
generating Nt electrons when starting from 100 pairs of electrons and ions. 

For the transfer from particle to hybrid simulation, we first determine the position of the 
model interface and the length of the buffer region, which is chosen dependent on the field ahead 
of the ionization front as we will discuss later Suppose the model interface is chosen at Zd - d*lSz 
and the buffer region is on the interval [zc, Zd\ as shown in Fig. [TO] where c,d e {1, 2, ■ • ■ , M.} 
and c < d. The particles in the future fluid region in the interval [zq, Zd] are then averaged in this 
particular time step to the densities on the grid G, while the particles outside the particle region 
and buffer region on the interval [zq, Zc] are removed from the particle list. Note that in the buff'er 
region [zc, Zd], the electron and ion densities are calculated on the underlying continuum grid 
while the spatial and kinetic information of the individual particles is also maintained. 

Now a general time step in the further evolution is described. At the beginning of a time step 
lasting from f„ to f„+i, we have the electric field E„, the electron and ion densities rig, rip in the 
fluid region and the positions of the particles in the particle and buffer region at time t„, and the 
velocities of these particles at time f„+i/2 (since with the leap-frog algorithm of the particle model, 
in the beginning of the time step, the positions x are known at t„ and the velocities v at tn+xji)- 
In step 1 (as shown in Fig.fTOli. the particles in the particle region and the buffer region are first 
moved to taking the stochastic collisions into account through the Monte Carlo algorithm. 
While updating the particle positions, the number of electrons crossing the model interface is 
recorded. 

Once all particle positions are updated to the next time step, in step 2 we also update the 
densities in the fluid model. The continuum equations for the particle densities are discretized 
with a finite volume method, based on mass balances for all cells in the fluid region on the grid 
G. Eq. O) is used to compute the density fluxes on the face of each cell in the fluid region, except 
the one at the model interface Zd where the flux from the particle model is used. The ionization 
rate is calculated with Eq. (O in the extended fluid model. Given the fluxes on the cell faces and 
the ionization rate, particle densities at each cell in the fluid region can be updated to the next 
step using Eq. ([T]i and Eq. (|2]i with the nonlocal parameter functions yu,,, D„ and a„. 

The densities are calculated using the third order upwind-biased advection scheme. The two- 
stage Runge-Kutta method could be used with the given particle flux (which is an average flux 
over [f„, f„+i]). This will allow computations with larger time step At, but it also requires two 
updates of the fluid densities and electric fields per time step. Since updating the fluid densities 
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needs the particle flux across the model interface as an input for the boundary condition, the 
two-stage Runge-Kutta method is replaced by the simpler forward Euler method, as long as 
the chosen time step is small enough. The time step and grid size in the hybrid simulation are 
chosen as the same as in the particle simulation and the fluid simulation, see Section 12.2.21 i.e., 
as Az = 2.3yL(m and At - 0.3ps, therefore the simulation results are compared with each other 
with the same numerical discretization errors. 

Now both the particle positions and the densities are known in the respective regions at time 
f„+i. In step 3, the densities in the particle region [zd, Zm-] are obtained by mapping particles to 
the grid G. How this is done, is discussed further below. The particle densities are then known 
everywhere on the grid G and the electric field at time f„+i can be calculated. 

In step 4, the position of the model interface is determined either from the particle density or 
from the electric field criterion. While the density profile and the electric field are updated from 
tn to tn+i, the model interface may have to move one cell forward or to stay still. With E'^ = -50 
to -200 kV/cm, the ionization front needs about 30 to 3 time steps to cross one cell. That is, the 
model interface will on average stay at one cell face for 3 to 30 time steps before it moves to the 
next cell face. If it stays still, the particles which fly out of the particle and buffer regions are 
removed from the particle list. If the model interface moves one cell forward, the fluid region 
is extended one cell into the particle region. Meanwhile, the buffer region also moves one cell 
forward and the region from Zc to becomes part of the fluid region, and particles there will be 
removed from the particle list. 

Finally we use E~ at f„+i to update the electron velocities inside the particle and buffer region 
to fn+3/2- This finishes one hybrid time step. 

We would like to remark here that one should be careful about the technique of mapping 
the particles to the densities in the particle region. One can use zero-order weighting by simply 
counting the number of particles within one cell, or first-order weighting (also called Particle In 
Cell or PIC) 1I54I1, which linearly interpolates charges to the neighboring cells, or some higher- 
order weighting techniques like quadratic or cubic splines. The particle model has been tested 
with various sizes of time steps and cell sizes. The numerical discretization errors converge 
to zero as time step and cell size decrease when using either the zero-order weighting or the 
first-order weighting, but the convergence was faster with first-order weighting. The first-order 
mapping is the most used technique in particle simulations for plasmas because: "As a cloud 
moves through the grid, the first-order weighting contributes to density much more smoothly 
than zero-order weighting; hence, the resultant plasma density and field will have much less 
noise and be more acceptable for most plasma simulation problems." js^l- However, in the 
hybrid computation, the zero-order weighting guarantees total charge conservation in the system, 
while the first-order weighting may cause charge loss or gain near the model interface when it is 
applied in a non-uniform density region. Therefore, zero-order weighting is implemented in our 
hybrid model. 

3.3. The position of the model interface 



We have discussed the suitable position of the model interface in 112811 where the particle 
model and the classical fluid model were coupled in space. It was shown that in the front part of 
the ionization front, the mean electron energies in the classical fluid model were lower than in the 
particle model, while behind the front the electron energies in both models are in good agreement. 
Therefore for a good agreement between hybrid and particle simulation, in the hybrid model the 
fluid model can be used behind the front, and a sufficiently large part of the density decay region 
should be covered by the particle model. 
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(a) Old hybrid model (b) New hybrid model 



Figure 11: The hybrid model (dotted) coupled with the classical fluid model (left column) or the extended fluid model 
(right column) is compared with the classical fluid model (dashed) and the particle model (solid) for a front propagating 
into a field of E'^ = -100 kV/cm. The model interface is located at three different levels of the field E = yE^ with 
y = 0.6, 0.9, and 0.98 shown in the upper, middle and lower panel, respectively. 
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Now within the hybrid model, the classical fluid model is replaced by the extended fluid 
model. Since an extra non-local term was introduced, such that the fluid simulation results are 
now closer to the particle results, we expect that a smaller region needs to be covered by the 
particle model without changing the result of the hybrid model. 

Fig. [TT] shows the simulation results of an ionization front propagating into a field of = 
-100 kV/cm; the left column shows results of the old hybrid model and the right column those 
of the new hybrid model. The model interface is located at three different positions: E - y 
where y - 0.6, 0.9, and 0.98, corresponding to the upper, middle and lower panel, respectively. 
Both in the old and the new hybrid simulation, a long buff'er region of 32 Az has been used to 
ensure the stable interaction of the models at the model interface. On the left, the old hybrid 
model for y - 0.6 generates the same electron density behind the front as the particle model, 
for y - 0.98 the same density as in the classical fluid model, and for y - 0.9 some intermediate 
density; we recall that the density difference between particle and classical fluid model is 20 
% for this field. On the right where the extended fluid model is used, the new hybrid model 
produces similar results as the particle model in all three cases y - 0.6, y = 0.9 and y - 0.98. 

The comparison shows that when the extended fluid model is used instead of the classical 
fluid model, the performance of the hybrid model is largely improved. Using the extended fluid 
model in the hybrid computation, the particle model can focus on a smaller portion of the front 
where the electron density is much lower while the electron density behind the front is still 
calculated with high accuracy. This greatly reduces the number of electrons that need to be 
followed in the particle region. It will give a substantial improvement of computational efficiency 
in a 3D simulation where millions of electrons will be pushed into the fluid region when the fluid 
model can be applied further ahead within the front. 

To determine the proper position of the model interface, different positions have been tested 
in planar front simulations. Fig. [T2] shows the relative density differences defined in Eq. ( fT4] i 
between the hybrid simulation and the particle simulation as a function of the position of the 
model interface; here the field ahead of the ionization front is E^ - -100 kV/cm. The model 
interface is placed at an electron density level tig = x iie^max where x varies. For x - \, the model 
interface is at the density peak ngj,,ax of the ionization front. From x - 1.0 to 0.1, the model 
interface moves forward within the front, and x - corresponds to a full fluid simulation. The 
full particle simulation is denoted as x = 1.1. 

3.4. The buffer region 

The particle model extends backward, beyond the model interface, into the buffer region 
where particle and fluid model coexist; it supplies particle fluxes for the fluid model at the model 
interface, as illustrated in Fig[TOl However, correct particle fluxes require correct particle statis- 
tics within the buffer region whose length should be as small as possible to reduce computation 
costs, but larger than the electron energy relaxation length |@]. 

3.4.1. Adding particles in the buffer region 

As we have mentioned, to ensure a stable electron flux at the model interface, in general new 
particles need to be introduced into the buffer region, that have to be drawn from appropriate 
distributions in configuration space. This would pose a particular problem for the streamer sim- 
ulation, since Maxwellian or even Druyvesteyn ISSf distributions are inaccurate. However, even 
for negative streamers, the electrons on average move somewhat slower than the whole ioniza- 
tion front, which means that the electrons on average are moving from the particle region into 
the fluid region. 
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Figure 12: The relative density differences Eq. U4t between the hybrid simulation and the particle simulation as a 
function of the position of the model interface = ,v for £^ = -100 kV/cm, The value ,v = accounts for the 

extended fluid model. The full particle simulation is denoted as x = 1.1. 



This is because in a pulled negative ionization front, the front velocity is determined in the 
leading edge of the front where the electric field and therefore the electron drift velocity is high- 
est. But the front velocity is determined by this maximal electron velocity plus an additional 
positive contribution of electron diffusion and impact ionization ilSII, therefore the front every- 
where is faster than the local mean electron velocity. 

Suppose that when the computation is changing from the pure particle regime to the hybrid 
regime, we create a buffer region which is long enough to relax most of the fast electrons. The 
particles in the buffer region are the heritage of the pure particle simulation and will keep being 
followed by the particle model. Across one end of the buffer region, the model interface at Zd, 
particles move freely. But across the other end of the buffer region at Zc, there are electrons flying 
out of the buffer cells but no electrons enter, i.e., when the model interface stays still, the buffer 
cells are losing particles. If we add electrons in the buffer region near Zc to create a flux from the 
fluid cell Zc-i into the buffer region Zc, although the particle loss is compensated, these inflowing 
particles have hardly any influence on the flux across the model interface when Zd-Zc is sufficient 
long. Once the model interface moves forward, those added particles near Zc will fall into the 
fluid region and be removed again. 

This hypothesis has been tested in hybrid simulations with two kind of fluxes used at Zc-\i2'- 

• No influx: electrons can only fly out over Zc-xji, but not fly back. 

• Reflected influx: electrons that fly out over Zc-i/2 Ay back immediately with inverted ve- 
locity. 

The hybrid simulation has been carried out for a planar front at -100 kV/cm with the model 
interface at £ = 0.6, 0.8, 0.9 E^. With a long buffer region of 10 Az, there is no influence of the 
inflowing particles. We even tested the artificial case of "double reflection influx", where twice 
as many electrons fly back with doubled energies. Even then, there was no notable influence. 

So we conclude that if electrons move on average more slowly than the front, the electron 
loss at the end of a sufficiently long buffer region does not affect the calculation of particle fluxes 
at the model interface. Therefore the particles lost at the end of the buffer region can be ignored 
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Figure 13: Hybrid calculation of an ionization front propagating into the field of —100 kV/cm with the model interface 
placed at 0.9 E*. The plots show the average electron flux (first panel) across the model interface, and densities (second 
panel) and mean energies (third panel) of the electrons at the model interface as a function of the buffer length region. It 
shows that all the quantities converge when buffer length increases. 



and new electrons do not need to be created artificially. This actually leads to a simpler technique 
of model coupling in which the correct energy distribution of particles is not a concern any more. 
This technique is not only suitable for streamer simulations, but for all two phase problems 
where the front speed is higher than the speed of the individual particles, as is generally the case 
in "pulled front" problems i33ll34ll . cf. section [T3] For example, bacterial growth and transport, 
or wound healing as a front propagation problem can be treated in this way. 



3.4.2. The length of the buffer region 

A reasonable particle flux can be obtained without artificially adding new electrons in a 
sufficiently long buffer region. But how long is enough long? Although electrons on average 
propagate slower than the model interface, the high energy electrons at the tail of the energy 
distribution may move faster than the ionization front in a short time interval. The buffer region 
should be long enough that most of the fast electrons can relax to mean energies within this short 
time. The relaxation length of the fast electrons at the ionization front has been discussed in |0]. 
It shows that three electron swarms with identical energies of 0.5, 5, or 50 eV, respectively, equi- 
librate within 2 ps through losing or gaining energy in collisions. If the buffer region is very 
short, high energy electrons might enter from the fluid region and contribute to the flux on the 
model interface, but since they are not be included in the buffer region, such short buffer region 
will result in a lower flux on the model interface. 

A stable flux can be obtained only when the buffer region is long enough to relax most high 
energy electrons to the local energies within the buffer region. An upper bound for the mean 
forward velocity (v;) of local electrons is the front velocity v/ * 7x10^ m/s at-lOOkV/cm. When 
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Table 2: The proper position of the model interfaces for different fields, where = x rig . 



an electron with 50 eV and with very small radial velocity components, i.e., with » v,-, relaxes 
to the local electron energy, its velocity in the forward direction v' decreases from ~ 4 x 10^ m/s 
to Vf within 2 ps, and it propagates 2 x 10"'^ ■ (v', + v/)/2 - 4.7 fim, which is approximately 2 
cells. 

Experiments with different lengths of the buffer region are shown in Fig. [13] Here we plot 
average flux, density and mean energy of the electrons at the model interface for buffer regions 
of different lengths. The field ahead the front is E'^ - -100 kV/cm and the position of the model 
interface is at £ = 0.9£'^. It is clear that as the length of the buffer region increases, not only the 
flux densities, but also the electron density and the mean electron energy converge to their limit 
values. Our computational cells are 2.3jt/m long, and a buffer region with the length of 2 cells 
can akeady give a reasonable flux at the model interface. 

4. Simulation results in different fields 

Having presented our hybrid simulations of the front propagating into a field of E^ = -100 
kV/cm in detail, we now summarize results for fields ranging from -50 to -200 kV/cm. 

In Fig. [141 we present the relative density discrepancies {ii^p„„ - In'g pan of the satu- 

rated electron density behind the ionization front for different fields; here hybr denotes hybrid 
simulations and part particle simulations. For each field, the model interface has been placed 
at electron density levels tie - x ne,max, with x - \ .0, 0.9, ...,0.1. A long buffer region with 
32 cells was used in the hybrid calculation to ensure a stable flux at the model interface. Two 
horizontal dashed lines have been added at +5% to assist the choice of the proper position of the 
model interface. For example for a field of E^ — -125 kV/cm, the relative density difference 
is limited to +2% even if only a very small part of the front x = 0.1 is covered by the particle 
model. But due to the flux convergence problem discussed in Sect. 13.41 the model interface will 
be put slightly more backwards to have a stable flux in our actual hybrid computational model, 
as the computational costs are actually determined by the position of the back end of the buffer 
region. 

We summarize the proper positions of the model interfaces for different fields in Table. |2l 
With the positions chosen from this table, we have tested different lengths of buffer regions. For 
the buffer lengths given, the electron flux, density and mean energy lie not more than 2% off the 
limit of a very long buffer region. 

5. Conclusion 

The particle model contains all essential microscopic physical mechanisms of a streamer 
ionization front, but the computational effort increases with the growing number of particles 
and rapidly exceeds computer memory. Gathering many real particles into one super-particle is 
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Figure 14: The relative density differences between the hybrid simulation and the particle simulation as a function of the 
position of the model interface in a range of fields. The model interface is at the position where the electron density level 
is 111 = xn^jnax- 

not an option, as it causes wrong statistics and causes numerical artifacts. The fluid model is 
much more eflicient than the particle model since the particles are approximated as continuous 
densities, but particle density fluctuations and large deviations of individual electron energies 
from the mean are not modeled. To combine the computational efficiency of the fluid model and 
the full physics of the particle model, a spatially hybrid model has been implemented for planar 
fronts. The hybrid model uses a particle model in the leading part of the ionization front where 
the total electron number is low and the field is high, and an extended fluid model for the rest. 
Once the scheme is extended to 3D, it can represent particle number fluctuations in the leading 
edge that might play a role in accelerating the branching process, or it can follow the run-away 
of individual electrons from the front as a possible cause of X-ray bursts and gamma-ray flashes. 

For constructing the hybrid model, the fluid model has been extended with a density gradient 
term in the impact ionization term. This extended fluid model has been introduced as an alter- 
native to the classical fluid model, since it maintains an electron flux that is consistent with the 
particle model. Both avalanche experiments and planar front experiments confirm that the ex- 
tended fluid model should be used in the hybrid calculation rather than the classical fluid model. 
(The fluid model with energy equation is discussed in appendix|A]as a possible alternative.) 

The hybrid algorithms and our numerical implementation are discussed in detail, in partic- 
ular, the position of the interface between particle and fluid model, and the construction of the 
buffer region at the interface. After testing the hybrid model at -100 kV/cm, we also present the 
hybrid simulation results for other fields, resulting in recommendations for the position of the 
model interface and the length of the buffer region in future 3D simulations. 

Although the spatially hybrid model is discussed in the context of streamer simulations, we 
would like to emphasize two particular points that are generic for pulled front problems. 1. 
When particle avalanches are used to calculate the transport and reaction coefficients for the 
fluid model, particle or fluid fronts agree particularly well, as they have the same front shape 
as a particle avalanche; this property holds only for pulled fronts and not for pushed or bistable 
fronts. 2. As the electrons on average move backward in a frame moving with the front velocity, 
the problem of creating particles with proper statistics in the buffer region does not need to be 
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solved. This problem here would be particular severe, as the distribution of electron energies 
and velocities are far from thermal and not even in equilibrium with the local electric field. But 
in many other cases, the distribution of inflowing particles is unavailable or difficult to obtain 
as well. The simplified coupling approach here off'ers an alternative for those problems without 
having to consider the influx on the back end of the bufifer region. 
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A. The electron energy equation 



A set of moment equations can be derived from the Boltzmann equation, in which the most 
significant quantities are the density, momentum, and energy of the electrons. The Bohzmann 
equation is of the form 



dt m \ . , 

where /(r, v, t) is the distribution function of the considered charged particle depending on time 
t, position r and velocity v, m is the particle mass, F is the electrical force exerted on the particle, 
and the term in the right hand side of Eq. (fT6b is the time rate of change of / due to collisions. 

The evolution of electron density, electron average momentum and electron mean energy are 
defined by the zero-, first-, and second-order moments of the distribution function li56Li57ll . For 
example, the density equations in the classical fluid model, Eq. ([TJ and Eq. (|2]i, can be deduced 
from the Boltzmann equation by integrating Eq. (fTST l over d\. The momentum equation can be 
obtained by multiplying Eq. (fTSI l by m\ and integrating over d\; The energy equation can be 
obtained by multiplying Eq. ( fT6] l by j v v and integrating over d\. 

Here we will mainly discuss the electron energy equation as a potential solution for the non- 
local effects reported in this paper and in [91]. We have shown in ^ that in the particle simulation, 
the mean energies of local electrons at the front are higher than in the classical fluid simulation, 
both in the electron avalanche and in the planar front, see figures 7 and 10 in jst]. In Sec- 
tion l2.2.2l we concluded that the local field and the local density approximation are not sufficient 
to describe the non-local spatial energy distribution appearing in the particle simulation, which 
causes the flux discrepancy between the particle and fluid description in a hybrid calculation, 
see Section IZ2] This problem is addressed by an extended fluid model which approximates the 
non-equilibrium ionization rate with a density gradient term, see Section l23] However, since the 
electron behavior (drift, diffusion, collision frequency) is determined by the electron energies, 
that are approximated by the local electric field, the equation for the local mean electron energy 
is here explored as an alternative to the gradient expansion. The energy equation can describe the 
non-equilibrium behavior of electrons in an electric field if ;) the transport and reaction coeffi- 
cients can be approximated as functions of the mean electron energies, and //) the electron mean 
energies are calculated with an electron energy equation. 

We first review the transport coefficients and redefine the electron flux jg and ionization 
source term S in the density equations ([T]i and (|2]i 

j, = -Me) E - D(e) ■ V«„ (17) 
S = KMe)E|a(e), (18) 

where the transport coefficients ju and D, and the ionization rate a are now functions of the mean 
electron energy e. 

The electron energy equation has been used in various plasma simulations, for example, for 
rf (radio frequence) discharges [Isil .SS. 59, 60J, for lamp breakdown processes [!6lll . and also for 



streamers |l62l|45l|52t]. The equation has the general form: 



+ V-J, = -ej,-E+ — — , (19) 



dt \ dt 

where is the energy flux (not to be confused with j^,), and (^^|j^)c represents the total energy 
loss due to collisions. 
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The energy flux can be written in drift-diffusion form similarly to the density flux as 



(20) 



where the fj.^ and are the energy rnobility and the energy diffusion coefficients. Many fluid 
models in the literature 1 62, 51, 58, 5^ 61, 631 use the energy transport coefficients given by 



and 



D, = - D. 

3 



(21) 



These approximations can be derived by assuming a Maxwellian EEDF, a constant momentum- 
transfer frequency and constant kinetic pressure ll52ll53ll . 

Concerning the energy collision term, the rate coefficients can be directly obtained from 
electron swarm experiments. The rate coefficients for excitational and ionization collisions are 
calculated in the form of Townsend coefficients ait{e), where k stands for the kth inelastical col- 
lision type with the energy loss Ck. If we have m electron-neutral inelastical collision processes, 
the time rate of change of electron energy due to collisions is the following: 



din^e) 
dt 



A practical form of the energy equation is, 

diriei) 5 

— + 3V.(j.e) 



(22) 



k=l 



-e - E - ^ Ijfl Ok ek. 



(23) 



k=i 



The energy equation coupled with the density equation is then solved for an electron avalanche 
The simulation results are presented in Fig. [15] As in Section l2.2.2l and Fig. HI a small electron 
swarm is first generated by the particle simulation at f = f i , then fluid model and particle model 
follow it to f = f2- Plotted are the same quantities as in Fig.|6]plus fluid and particle results for 
the local mean electron energies at time f2, the particle results for the electron energy were first 
reported in Fig. 10 in 1@1. 

As shown in Fig. [15] the fluid model with energy equation generates a similar swarm as the 
particle model while small differences remain in velocity and maximum electron density. The 
local energy description in particle simulation and fluid simulation agrees surprisingly well. Note 
that this nice agreement is obtained when the energy flux coefficients adopt the commonly used 
approximation as shown in Eq. dTTT i. where the electron energy distribution function is assumed to 
be Maxwellian, though previous studies have shown that the electron energy distribution function 
is not Maxwellian during an avalanche in the high field, cf. Fig. 3 in 1@1. 

The energy equation is an alternative for the gradient expansion or vice versa, since both 
methods can describe the deviation of electrons from the local field approximation. The deviation 
has been addressed as the 'non-local' effect in the density and field gradient approach l47ll . or as 
the 'non-equilibrium' effect in the multi-moment equation approach 15^ 45, 521, which shows 
the two different perspectives of these two approaches. 

Which approach is better in the hybrid model? The accuracy of the coefficients plays an es- 
sential role. The electron 'non-local' rate coefficients ki in Eq. (O, including all combinations of 
the density gradients and field gradients, have been calculated by solving the Boltzmann equation 
in a two-term approximation ll46ll . Some authors working with the energy equation i63ll64ll53[l 
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Figure 15: Data of the particle model and of the fluid model with energy equation are shown: the density profile of the 
electron swarms and the initial present electrons in a field of -100 kV/cm are presented at two dilferent time ti and t2', we 
also show the mean energies of electrons of swarms at time t2. 

also avoid the commonly used energy transport coefficients approximation in Eq. ( 1211 1. and cal- 
culate them more precisely using the two-term Boltzmann equation approximation. However, it 
is not clear, how consistent the calculated coefficients from the two-term Boltzmann equation are 
with the particle description, while consistency is important for the hybrid coupling. When the 
classical fluid model is extended with the first and most important term in Eq. (|7]l, the density 
gradient, the rate coefficient of this term can be easily calculated from the particle swarm experi- 
ment. Since this simple extension already brings a good agreement in electron flux and since it is 
sufficient for the hybrid coupling, more gradient terms or the energy equation are not introduced 
in the fluid model in this paper. 
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